/* 3D Navier-Stokes Solver -> Level Set -> Marching Cubes */

#ifndef NAVIER_H
#define NAVIER_H

#include <iostream>
#include <vector>
#include "typen.h"
#include "matrix.h"
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include "velocity.h"

//change the following definition if you have too much objects in your NAV file
typedef unsigned int flagtype;
const flagtype flagtypemsb=UINT_MAX & ~(UINT_MAX >> 1);

/* Class Navier contains all kind of lists needed for different boundary treatments for 
   the velocities, pressure, temperature and chemicals. */    
class Navier
{
   protected:
      Matrix <double>  P; // P is the array for pressure
      Matrix <double>  T; // T is the array for temperature
      Matrix <double>  BiCG_r0, BiCG_rj, BiCG_pj, BiCG_vj, BiCG_sj, BiCG_tj, MatBuf;// Temporary array needed for the BiCGSTAB solver scheme 
      Matrix <double>    KoeffPreCond;// KoeffPreCond array contains the coefficients for the jacobi preconditioner included in the BiCGStab scheme 
      unsigned save_mem;
      Matrix <double> RHS;// RHS is the right hand side array of the poisson equation      
      Matrix <double> F, G, H, T2, Ft, Fl;// some temporarly used arrays for velocity, temperature and level-set
      //NS_List<int> *listAirToFluid;
      std::vector <int> listAirToFluid;
      std::vector <int> listFluidToAir;
      std::vector <int> listAirBorder;
      //std::vector <int> slip_list; 
      //std::vector <int> outflow_list;
      int num_iter_Navier;
      int slip_list[10000];
      int outflow_list[10000];
      //int inflow_list[10000];
      int slip_list_size;
      int outflow_list_size;
      //int inflow_list_size;
      int iterations, alliterations;//number of iterations for the poisson solver
      int n; //number of steps
      int actual_fgh;// tag which defines if F,G,H arrays are from actual or previous time step
      bool flagsixfluidneighbors;
      void DefinirCenaInicial(int argc,char* argv[]) ;
      void Create_slip_list(int &count);
      void Create_outflow_list(int &count);
      void Prepare_Air_To_Fluid(int *DirectList,int *size_updated_List);
      void orderAirToFluid(int *list, int sizeList, unsigned int *flaglevel, int *DirectList,int *size_updated_List);
      void Initializer_No_Face_or_Edge_Neighbor( int i, int j, int k);//Esse Inicializador quando  ac�lula n�o tem vizinhos fluido
      void NewFluidCellInitializer(int i, int j, int k);//Inicializador das c�lulas que eram Air e passaram a ser fluido
      void Prepare_Structure_for_Slip_Outflow_Cells(/*std::vector <int>*/int *list, int &slip_list_size, int list_offset,
                                                    int list_two_free_faces);
      void Init_Velocity_Inflow();
      void Call_Complete_Cell();
      void Complete_Cell(int i,int j, int k, int freefacespar);
      void Call_Update_Slip_Outflow_Cells();
      void Update_Slip_Outflow_Cells(int i1, int j1, int k1, int i2, int j2, int k2,int *factor0,int *factor1);
   
      void AdapUVW(); //velocity correction with pressure gradient due to chorin projection method (see User's Guide)
      void CompFGH(Velocity gridCellUVW);//calculate temporary velocity fields (in space) due to chorin projection method (see User's Guide)
      
      void CompF(int i, int j, int k);
      void CompG(int i, int j, int k);
      void CompH(int i, int j, int k);
      void CompFM1(Velocity gridCellUVW, int i, int j, int k,int gama);
      void CompGM1(Velocity gridCellUVW, int i, int j, int k,int gama);
      void CompHM1(Velocity gridCellUVW, int i, int j, int k,int gama);
      
      double BiCGStab();//BiCGStab poisson solver method 
      //!calculate temporary velocity fields (in time) due to chorin projection method (see User's Guide)
      void CompTUVWfromFGH(int timestepmethod); // computes \tilde u calcula u til      
      void CompRHS();// compute right hand side of poisson equation
      double TimeStep();// calculate new time step
      //void CompTG(Matrix<double>& CH, double alpha, double diffconst);//compute transport equation in space
      double CompPoisson(/*double starttime*/);//compute poisson equation
      double Poisson();//compute poisson equation
      void CalcKoeffPreCond();//compute coefficients for jacobi precondidioner
      void SetPBorder(Matrix<double>& B,int cg);//set color dependent pressure boundary conditions (done in every SOR iteration)
      void SetPALLBorder(Matrix <double> &X);//set global pressure boundary conditions
      void SetPLocalBorder(Matrix <double> &X, int i, int j, int k);//set local pressure boundary conditions    
      void SetObstacleCond(int method);//set velocity/temperature/chemical boundary conditions
      void SetObstacleCondTilde(int method);
      //! set boundary conditions for transport equation
      void SetObstacleForTG(Matrix<double>& M, int n);
      double CalcRes();              
      //  the following functions are used locally
      //interpolate values by third order lagrange interpolation
      double LagrangeIntp(double LL, double L, double R, double RR, double DD, double D, double P, double PP);
      double Precondition(double P, int i, int j, int k);//compute the preconditioning 
      double PreconditionBack(double P, int i, int j, int k);//reverse computation of the preconditioning
      double MatVecMult(Matrix<double>& A,int i, int j, int k, int RES, bool flag_bool);//matrix vector multiplication
      //compute the mixed partial derivative
      double DUV(double LL,double L,double M,double R,double RR,double KL,double KR,double DD,double DM,double D,double DP,double PP,double alpha,int border);
      //compute the partial derivative
      double DUU(double LL,double L,double M,double R,double RR,double DD,double D,double DP,double PP,double alpha, int border);
      double DL(double LLL, double LL, double L, double M, double R, double RR, double RRR, 
                  double DDD, double DD, double D, double DP, double PP, double PPP, 
                  double sign, int reinit);  
          
   public://as defini como p�blico para poder atribu�las direatamente a gridU,gridV,gridW
      Scene S; // S contains all physical and computational parameters, the so called scene description physical+computational parameters
      Matrix <double>  U; // U is the array of the velocity component in x-direction
      Matrix <double>  V; // V is the array of the velocity component in y-direction
      Matrix <double>  W; // W is the array of the velocity component in z-direction
      Matrix <double>  Phi; // Valor de phi para ajudar a verificar se uma determinada c�lula � fluido ou n�o 
      Matrix <flagtype>   flag; //flag field, this array contains the geometry and tags all cells types (FLUID, AIR,...)  
      //int *finalList;
      double UVWmax;
      void  DoIt(Velocity gridCellUVW);
      Navier(int argc,char* argv[]);
      ~Navier();
};

#endif
